/*
  Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.
  
*/

/// \file
/// Contains the grid_cell class.
// $Id$


// (This file is included automatically by grid.h)

#include "blitz/numinquire.h"
#include "ray.h"

/** A cell in a hierarchichal grid. The grid_cell objects are
    contained in a grid represented by a grid_base object. A cell
    either contains a subgrid or, if it is a leaf cell, a T_data
    pointer to whatever is being contained in the cells. */
template<typename cell_data_type> class mcrx::grid_cell {
public:
  typedef mcrx::T_float T_float;
  typedef grid_base<cell_data_type> T_grid;
  typedef grid_cell T_cell;
  typedef cell_data_type T_data;
  typedef unsigned int T_indextype;
  typedef blitz::TinyVector<T_indextype, 3> T_index;
private:
  typedef int T_neighbor_offset;
  friend class blitz::TinyVector < T_cell, 8 >;
  friend class blitz::MemoryBlock<T_cell> ;
  friend class grid_base<cell_data_type> ;
  /// Copy constructor does NOT copy sub grids or cell data!
  grid_cell (const T_cell&);
  T_cell& operator= (const T_cell&);
protected:

  /// Bit set indicating whether the neighbor pointer is bad.
  mutable short int neighbor_set_bad;
  bool leaf; ///< Indicates whether cell is refined or not.

  /// Index vector of this cell in its grid
  T_index index;

  union {
    /// Pointer to a possible refined grid in this cell (mutually
    /// exclusive with celldata).
    T_grid* sub;
    /// Pointer to data for this cell (mutually exclusive with sub)
    T_data* celldata;
  };

  /// Pointer to the grid the cell is in.
  T_grid* host;

  /// \name Caching stuff for the neighbor-finding routines.
  ///@{
  /// Cache vector containing pointers to neighboring cells.
  mutable blitz::TinyVector <T_neighbor_offset, 6 > neighbors;
  const T_cell* slow_neighbor(const vec3d &,int, bool,int) const;
  const T_cell* neighbor_samelevel(const vec3d &,int, bool,int) const ;
  ///@}

  void placement_new_destructor ();

  // Neighbor search and ray/cell intersections
  const T_cell* neighbor (const vec3d& p,int dim, bool positive_dir) const;
  /// Finds neighboring cells. Just calls the const version.
  T_cell* neighbor (const vec3d&p ,int dim, bool positive_dir)  {
    return const_cast<T_cell*> 
      (const_cast <const T_cell*> (this)->neighbor(p, dim, positive_dir));};

public:
  grid_cell ();
  // creates a cell with the specified index in the specified host grid
  grid_cell(T_grid *hg, T_index indx);
  ~grid_cell();

  /// Returns the pointer to the cell data object.
  T_data* data() const {assert (leaf); assert (celldata); return celldata;};
  /// Sets the celldata pointer, deleting the object previously
  /// pointed to first.
  void set_data(T_data* d) {assert (leaf); if (celldata) delete celldata;
  assert (d); celldata = d;};
  /** Sets the celldata pointer, without deleting the object
      previously pointed to. This is used when the current celldata
      was created using placement (this is not optimal). */
  void set_data_nodelete(T_data* d) {assert (leaf); celldata = d;};

  class cell_tracker;

  // Finds exit point of a ray, assuming it is inside the cell.
  T_float
  intersection_from_within (const ray_base&, T_float, 
			    T_densities&, cell_tracker&);
  /** Finds the fractional distance along a ray that corresponds to a
      certain column density through the cell. Since the cells are
      uniform density, this is trivial for this grid. */
  T_float column_to_distance_fraction(/// The column density fraction you want distance for.
				      T_float fn,
				      /// The starting position
				      const vec3d& pos, 
				      /// The direction along which the column is calculated
				      const vec3d& dir,
				      /// Total length of the segment
				      T_float dl) const {
    return fn; };

  // refinement
  void refine ();
  void refine (void*&);
  void unrefine (T_data*const);
  const T_grid* host_grid () const {assert (host); return host;};
  T_grid* host_grid () {assert (host); return host;};
  const T_grid* sub_grid () const {assert (!leaf); assert (sub); return sub;};
  T_grid* sub_grid () {assert (!leaf); assert (sub); return sub;};

  // simple data retrieval
  /// Returns whether the cell is a leaf cell or not.
  bool is_leaf () const {return leaf;};
  /// Returns the index vector in the grid of this cell.
  const T_index& getindex () const {return index;};
  /// Returns the min point of the cell boundary.
  const vec3d getmin() const {
    return host_grid ()->getmin() + getindex ()*getsize();};
  /// Returns the max point of the cell boundary.
  const vec3d getmax() const {
    return getmin () + getsize ();};
  /// Returns the size of the cell
  const vec3d& getsize() const {return host_grid ()->getsize ();};
  /// Returns the center point of the cell
  const vec3d getcenter() const{
    return host_grid ()->getmin() + 
      (getindex ()+vec3d (.5, .5, .5))*getsize();};
  /// Checks whether a position is within the cell.
  bool contains(const vec3d& p) const {
    return all (p >= getmin ()) && all (p < getmax ());};
  /// Returns the volume of the cell.
  T_float volume() const {
    const T_float v=product (getsize ());assert(v>0); return v;};

  /// Returns the top-level grid of the cell.
  const T_grid* top_grid() const {
    if(host->parent) return host->parent->top_grid(); else return host;}
};


/** This class abstracts the concept of a "current cell" from the
    radiation transfer code. This is necessary because in general we
    can't assume a 1:1 correspondence between ray steps and grid cells
    (where a grid cell is something that has e.g. a density and
    intensity). In the case of this grid, it is just a trivial
    forwarder. */
template<typename T> class mcrx::grid_cell<T>::cell_tracker {
public:
  typedef grid_cell<T> T_cell;

private:
  T_cell* cell_;

public:
  cell_tracker() : cell_(0) {};
  cell_tracker(T_cell* c) : cell_(c) {};

  bool operator==(const cell_tracker& rhs) const {return cell_==rhs.cell_;};

  /// Set tracker to invalid. (Used to be just a simple assignment in xfer.)
  void reset () {cell_=0;};

  /// Nop for this class.
  void reset_intersections() const {};

  /// Converts to a bool indicating whether cell is valid or not.
  operator bool() const {return cell_ != 0;};

  T_cell* cell() const {return cell_;};

  T_float
  intersection_from_within (const ray_base& r, T_float max_len, 
			    T_densities& dn) {
    return cell_->intersection_from_within(r, max_len, dn, *this); };

  T_float column_to_distance_fraction(T_float fn,
				      const vec3d& pos, 
				      const vec3d& dir,
				      T_float dl) const {
    return cell_->column_to_distance_fraction(fn, pos, dir, dl); };

  bool contains(const vec3d& p) const {
    return cell_->contains(p); };

  void assert_contains(const vec3d& p) const {
    assert(contains(p)); };

  T_float volume() const {
      return cell_->volume(); };

  typename T_cell::T_data* data() const {
    return cell_->data(); };
};


// *** constructors/destructors ***

/// Default constructor is just for the array allocations.
template < typename cell_data_type >
mcrx::grid_cell<cell_data_type>::grid_cell():
  neighbor_set_bad (0), leaf (true), index (0,0,0),
  celldata (0), 
  host (0), neighbors (0)
{
  assert (sizeof (neighbor_set_bad) == 2);
  assert(sizeof(long int) >= sizeof(void*));
}

template < typename cell_data_type >
mcrx::grid_cell<cell_data_type>::grid_cell(T_grid *hg, T_index indx):
  neighbor_set_bad (0), leaf (true), index (indx), celldata (0), 
  host (hg), neighbors (0) 
{
  assert (sizeof (neighbor_set_bad) == 2);
  assert(sizeof(long int) >= sizeof(void*));
}

/** copy constructor ONLY copies host grid and index, not any possible
    sub grids or cell data.  These fields are initialized to 0.  This
    is unfortunate but we need the copy constructor for the grid
    constructors. */
template < typename cell_data_type >
mcrx::grid_cell<cell_data_type>::grid_cell(const T_cell& c):
  leaf (true), celldata (0), host (c.host), index (c.index),
  neighbors (0), neighbor_set_bad (0)
{
}

template < typename cell_data_type >
mcrx::grid_cell<cell_data_type> ::~grid_cell()
{
  if (is_leaf ()) {
    if (celldata)
      delete celldata;
  } 
  else
    if (sub) 
      delete sub;
}

/** "Destructor" when placement new is used. The idea here is that we
    want to preempt the actual destructor from deleting sub if we have
    used placement new, so the top-level (adaptive_grid) destructor
    checks for that and calls this function first if that is the
    case. */
template < typename cell_data_type >
void mcrx::grid_cell<cell_data_type>::placement_new_destructor ()
{
  if (!is_leaf ()) {
    assert (sub);
    sub->placement_new_destructor();
    sub->~grid_base<cell_data_type>();
    sub = 0;
  }
}

// *** grid_cell methods ***

template<typename cell_data_type>
mcrx::grid_cell<cell_data_type>&
mcrx::grid_cell<cell_data_type>::operator= (const T_cell& rhs) 
{
  if (this == &rhs)
    return *this;

  // assignment operator does NOT copy sub grids/cell data, so it
  // should never be used for this.  Basically it's only here because
  // the TinyVector constructor (in octogrid) needs it.  (This sucks
  // because it opens a gaping hole in the management of the sub
  // grid/cell data pointers.)
  assert (rhs.sub == 0);
  
  leaf = rhs.leaf;
  assert (!sub);
  sub = rhs.sub; assert(!rhs.sub);
  host = rhs.host;
  index = rhs.index;
  neighbors = 0;
  neighbor_set_bad = 0;

  return (*this);
}


/** Finds neighboring grid cell. This is done by projecting position p
    in dimension dim, direction dir (-1 lower end, +1 upper end, other
    results undefined!)  if p is not within this cell weird things may
    or may not happen, but the call doesn't make sense then
    anyway. (For non-refined grids p doesn't matter.) If a cached
    result is available, it is used, otherwise slow_neighbor is
    called. */
template < typename cell_data_type >
inline const typename mcrx::grid_cell<cell_data_type>::T_cell*
mcrx::grid_cell<cell_data_type>::neighbor (const vec3d& p,
                                           int dim, bool positive_dir) const
{
  assert ((dim >= 0) && (dim < 3));
  const T_float tol = 1e-10;
  
  // First check if we have a cached neighbor in this direction
  const int dirindex=(dim<<1) + (positive_dir?  1:0); 
  const T_cell* returnval;

  if( (neighbor_set_bad & (0x101 << dirindex)) == (0x001 << dirindex)) {
    // Awesome - we have done this before, so we can (almost) just return
    // the pointer.  But first some checks:

    // Extract the pointer from the offset
    const T_cell* const nc= 
      neighbors [dirindex] != 0? reinterpret_cast<T_cell*>
      (reinterpret_cast<long int> (this) + neighbors[dirindex]): 0;
    if( nc && !nc->is_leaf()) {
      // The cell exists, and is subdivided, 
      // we must call locate to find the actual leaf grid cell

      // check that the caching returned something sensible
      assert(positive_dir? 
	     (std::abs(getmax()[dim] - nc->getmin()[dim]) < tol) : 
	     (std::abs(getmin()[dim] - nc->getmax()[dim]) < tol));
      // and that the dimensions are unchanged in other dims
      assert((dim ==0) || ((nc->getmin()[0]-tol <= p[0]) &&
			   (nc->getmax()[0]+tol >= p[0])));
      assert((dim ==1) || ((nc->getmin()[1]-tol <= p[1]) && 
			   (nc->getmax()[1]+tol >= p[1])));
      assert((dim ==2) || ((nc->getmin()[2]-tol <= p[2]) &&
			   (nc->getmax()[2]+tol >= p[2])));

      const T_cell* const nc2 = 
	nc->sub_grid()->locate(p ).cell();
      returnval = nc2;
    }
    else {
      // It's not subdivided, or simply doesn't exist - just return it
      returnval = nc;
    }
  }
  else
    // we don't have a cached neighbor pointer, do it the slow way
    returnval = slow_neighbor (p,dim, positive_dir,dirindex);

  // To make dead sure this routine is working:
  // check that cell is actually bordering in dim
  assert(positive_dir? 
	 !returnval || 
	 (std::abs(getmax()[dim] - returnval->getmin()[dim]) < tol) :
	 !returnval || 
	 (std::abs(getmin()[dim] - returnval->getmax()[dim]) < tol));
  // and that the dimensions are unchanged in other dims
  assert((dim ==0) || !returnval || ((returnval->getmin()[0]-tol <= p[0]) &&
				     (returnval->getmax()[0]+tol >= p[0])));
  assert((dim ==1) || !returnval || ((returnval->getmin()[1]-tol <= p[1]) && 
				     (returnval->getmax()[1]+tol >= p[1])));
  assert((dim ==2) || !returnval || ((returnval->getmin()[2]-tol <= p[2]) &&
				     (returnval->getmax()[2]+tol >= p[2])));
  // and that it has a data object, if it's a valid cell
  assert((!returnval) || returnval->data());

  return returnval;
}

/** Finds neighboring grid cell by looking through the hierarchy of
    grids. Updates the vector of cached neighbors with the results so
    it's not called again. */
template < typename cell_data_type >
const typename mcrx::grid_cell<cell_data_type>::T_cell*
mcrx::grid_cell<cell_data_type>::slow_neighbor (const vec3d& p, int dim, 
                                                bool positive_dir, int
                                                dirindex)
const
{

  // We did not have a cached neighbor for that direction - now we must
  // perform the normal (slow) neighbor search and cache the result for 
  // future use

  // Get indices of this grid cell
  T_index ind = getindex ();

  // calculate the indices of the neighboring cell
  ind[dim]+= positive_dir?  1:-1;

  // check whether indices is within range
  const int nofdim= host_grid ()->getn() [dim];
  const int indofdim=ind[dim];
  const T_cell* nc;
  
  // now we have 3 possibilities: normal grid cell, refined grid cell or
  // outside the grid (in which case it may mean outside grid volume or
  // in another higher-level cell
  if( (indofdim<0) || (indofdim>=nofdim) )
    // outside this grid
    if(!host_grid()->parent_cell()) {
      // The grid we're in is the top-level grid, so the point must be
      // outside the entire grid volume - return null
      nc = 0;
    }
    else 
      // there is a higher-level grid - find the neighboring cell in that grid
      // Recursively call neighbor_samelevel for the higher-level cell until 
      // we return with the neighboring cell of the same refinement level
      //nc= host_grid ()->parent_cell()->neighbor_samelevel (p,dim,positive_dir,1);

      // to be sure we don't round to a diagonally opposed cell, we
      // use the cell center as coordinates for the neighbor search
      // here (we know the neighbor at the same level will have the
      // same size, so we don't need the exact position)
      nc= host_grid ()->parent_cell()->neighbor_samelevel (getcenter(),
							   dim,
							   positive_dir,1);

  else 
    // So it is a valid grid cell in this grid
    // Get a pointer to it
    nc = host_grid ()->get_cell(ind );

  // (Do we need a mutex here?  I don't think race conditions matter
  // since two competing threads will try to write exactly the
  // same contents. The important thing is that neighbor_set will
  // never be true unless neighbors has already been filled, and
  // once true will never change.)
  bool bad = false;
  if (nc) {
    const long int offset = reinterpret_cast<long int> (nc) -
      reinterpret_cast<long int> (this);
    // Does the offset fit?
    if ((offset < blitz::huge (T_neighbor_offset ())) && 
        (offset > blitz::neghuge (T_neighbor_offset ())))
      // yes, store 
      neighbors [dirindex]= offset;
    else
      // no, mark as bad
      bad = true;
  }
  else {
    neighbors [dirindex] = 0;
  }

  neighbor_set_bad |= ((1 | (bad?0x100:0)) << dirindex);


  // now check that this actually worked

  // check that bad was set properly
  assert(bad == ((neighbor_set_bad & (0x100 << dirindex)) != 0));
  // check that set was set
  assert(neighbor_set_bad & (1 << dirindex));
  // and finally check that we get the right cell back
  assert( !((neighbor_set_bad & (0x101 << dirindex)) ==
	    (0x001 << dirindex)) || 
	  ((neighbors [dirindex] != 0? 
	    reinterpret_cast<T_cell*> (reinterpret_cast<long int> (this) + 
				       neighbors[dirindex]) :
	    0) == nc));

  // Now this becomes a possible race condition!  The access to 
  // neighbor_set is no longer atomic, so conceivably (but highly
  // unlikely) two threads can try to write at the same time.  If they
  // are writing the same bit, nothing will happen, but if they are
  // writing different bits only one of them will end up getting set.
  // However, this will only lead to another call to slow_neighbor,
  // nothing bad will happen.  However, it IS important that both the
  // set and bad bit are set at the same time, because otherwise
  // things can get fucked up.

  // But that cell may be furtherly refined as well, so to fix that
  // we do exactly what standard neighbor () does here
  if (nc &&!nc->is_leaf())
    return nc->sub_grid()->locate(p ).cell();
  else
    return nc;
}  

/** Finds the neighboring cell of the same refinement level as this
    cell. If the refinement subdivisor is different in the two cells
    we have no idea what will happen when we try to use it...  If the
    neighbor is less subdivided we just return the leaf cell. */
template < typename cell_data_type >
const typename mcrx::grid_cell<cell_data_type>::T_cell*
mcrx::grid_cell<cell_data_type>::neighbor_samelevel(const vec3d& p, int dim,
                                                    bool positive_dir,int
                                                    level)
  const
{
  // Get indices of this grid cell
  T_index ind = getindex ();

  // calculate the indices of the neighboring cell
  ind[dim]+= positive_dir?  1:-1;

  const int nofdim=host_grid ()->getn() [dim];
  const int indofdim=ind[dim];

  // now we have 3 possibilities: normal grid cell, refined grid cell or
  // outside the grid (in which case it may mean outside grid volume or
  // in another higher-level cell
  if( (indofdim<0) || (indofdim>=nofdim) )
    // outside this grid
    if(host_grid ()->parent_cell()) 
      // there is a higher-level grid - find the neighboring cell in that grid
      // Recursively call neighbor for the higher-level cell
      return host_grid ()->parent_cell()->
        neighbor_samelevel(p,dim, positive_dir,level+1);
    else
      // The grid we're in is the top-level grid, so the point must be
      // outside the entire grid volume - return null
      return 0;
  else {
    // So it is a valid grid cell in this grid
    // Get a pointer to it
    const T_cell* const nb = host_grid ()->get_cell(ind);
    
    // If it's not refined it's easy, just return it.
    if(nb->is_leaf()) 
      return nb;
    else {
      // It IS refined, so we call locate in that subgrid
      const T_cell*const nc = nb->sub_grid ()->locate_samelevel(p,level-1);
      // we should never ever get 0 back here because that indicates a
      // failure of locate
      assert (nc);
      return nc;
    }
  } //else
}  


/** Refines a grid_cell to contain a sub-grid. */
template < typename cell_data_type >
void mcrx::grid_cell<cell_data_type>::refine () 
{
  assert (is_leaf ());
  if ( is_leaf ()) {
    if (celldata)
      delete celldata;

    sub = host_grid ()->create_sub(this);
    leaf = false;
  }
  // should we do something drastic if we try to refine an already
  // refined cell?
}

/** "Placement new" version of refine. Refines a grid_cell and puts
    the sub grid in the specified location in memory. */
template < typename cell_data_type >
void mcrx::grid_cell<cell_data_type>::refine (void*& placement) 
{
  assert (is_leaf ());
  if ( is_leaf ()) {
    if (celldata)
      delete celldata;

    sub = host_grid ()->create_sub(this, placement);
    leaf = false;
  }
  // should we do something drastic if we try to refine an already
  // refined cell?
}

/** "Unrefines" a cell and replaces it with the cell data object d. */
template < typename cell_data_type >
void mcrx::grid_cell<cell_data_type>::unrefine (T_data* const d) 
{ 
  if (!is_leaf ()) {
    delete sub_grid ();
    leaf = true;
    celldata = d;
  }
  // what about doing something drastic if we try to unrefine a leaf cell?
}


/** Finds the distance the ray will travel before exiting the cell.
    The function assumes that the ray is currently inside the cell, it
    won't work correctly for other cases.  */
template < typename cell_data_type >
mcrx::T_float
mcrx::grid_cell<cell_data_type>::
intersection_from_within (const ray_base& ray, T_float max_len, 
			  T_densities& dn, cell_tracker& next)
{
  // Find the cell boundary in the direction of travel
  const vec3d curmin=getmin();
  const vec3d curmax= curmin+getsize ();
  const vec3d idir = ray.inverse_direction();
  // it's better here to compare the sign of *idir* instead of dir,
  // because dir can be -0 which will be treated as +0 and lead to an
  // incorrect sign in dl3 below. With idir those are +-inf which will
  // compare correctly for sign.
  const vec3d bound((idir [0]<0)?curmin[0]:curmax[0],
                    (idir [1]<0)?curmin[1]:curmax[1],
                    (idir [2]<0)?curmin[2]:curmax[2]);

  // Compute the line element to the grid boundary. Get the distance
  // to the cell boundaries. Here we can get NaN if the ray is exactly
  // on the boundary and parallel to it, because idir will be Inf and
  // bound-pos will be zero. In that case we set it to Inf indicating
  // that it should never cross the boundary.
  const vec3d dl3=where(abs(idir)<blitz::huge(T_float()),
			(bound -ray.position())*idir,
			blitz::huge(T_float()));
  assert(all(dl3 > -1e-10 || dl3!=dl3)); 

  // Find the *closest* (nonnegative) boundary
  // (negative numbers can happen due to rounding, but they should be small)
  // b contains the dim# of the closest boundary, dl the distance. 
  //const int b = (dl3[0]<dl3[1])?  (dl3[0]<dl3[2]?0:2): (dl3[1]<dl3[2]?1:2);
  const int b = ((dl3[0]<dl3[1]) && (dl3[0]>=0)) ?
    ( ((dl3[0]<dl3[2]) && (dl3[0]>=0)) ? 0:2) :
    ( ((dl3[1]<dl3[2]) && (dl3[1]>=0)) ?1:2);

  // Get distance from current position to the next intersection point
  // (or to the max length point if we will pass it)
  const T_float dl = std::min(max_len - ray.length(), 
			      dl3[b]);

  dn = dl*data()->get_absorber().densities();

  if (dl == max_len - ray.length()) {
    // hit max length, invalidate tracker
    next = 0;
  }
  else {
    next = neighbor(ray.position()+ray.direction()*dl3[b], b, idir[b]>0);
  }

  
  return dl;
}
